function quadrantMaximumNorm = positionCP(images, numImgToProcess, cpNum)
% cpNum           - number of control points
% numImgToProcess - number of images to process

%get the mean image 
meanI = getMeanImage(images);

[h w d] = size(images);

optimCPResults    = zeros(2, cpNum*cpNum, numImgToProcess);
optimAlphaResults = zeros(2, cpNum*cpNum, numImgToProcess);

% figure(1);
%imagesc(meanI); colormap gray;
%hold on;
%loop over the images
for i = 1:numImgToProcess
    source = meanI;
    target = images(:,:,i);  
    
    [deformedSrc optimCP optimAlpha] = pairRegistration(source, target);

    optimCPResults(:,:,i)    = optimCP;
    optimAlphaResults(:,:,i) = optimAlpha;
    
    %plot control points
    %quiver(optimCP(1,:), optimCP(2,:), optimAlpha(1,:), optimAlpha(2,:));
    %pause(0.01);
end
% hold off;

%loop through every control point
% figure(2);
% imagesc(meanI); colormap gray;
% hold on;
% for i = 1:cpNum*cpNum
%     max    = [NaN, NaN, NaN, NaN]; 
%     maxInd = [NaN, NaN, NaN, NaN]; 
%     
%     %loop through every image match
%     for j = 1:numImgToProcess
%         % vector V = <v,u> is from A to B, 
%         % V = < b1 - a1 , b2 - a2 >, ||V|| = sqrt(v1^2 + v2^2), angle = tan-1(v/u) in degrees
%         
%         A = optimCPResults(:,i,j);
%         B = optimAlphaResults(:,i,j);
%         V = B;       
%         
%         magV  = sqrt(V(1)^2 + V(2)^2);
%         angle = radtodeg(atan2(V(2),V(1)));
%                 
%         %divide into quadrants
%         if angle < 90 && angle >= 0
%             q = 1;
%         elseif angle < 180 && angle >= 90
%             q = 2;
%         elseif angle < -90 && angle >= -180
%             q = 3;
%         elseif angle < 0 && angle >= -90
%             q = 4;        
%         end
%         
%         %get the maximum in each quadrant
%         if isnan(max(1,q)) || max(1,q) < magV 
%             max(1,q)    = magV;
%             maxInd(1,q) = j;
%         end
%     end
%         
%     %plot max vector in each quadrant
%     total = 0; 
%     for q = 1:4
%         ind = maxInd(1,q);
%         
%         if isnan(ind)
%             continue;
%         end        
%         cp    = optimCPResults(:, i, ind);
%         alpha = optimAlphaResults(:, i, ind);
%         total = total + norm(alpha);
% %         quiver(cp(1,:), cp(2,:), alpha(1,:), alpha(2,:));        
% %         pause(0.01);
%     end    
% end
% hold off

quadrantMaximum = zeros(size(meanI));
quadrantMaximum = std(optimAlphaResults,1,3);
quadrantMaximumNorm = sqrt(quadrantMaximum(1,:).^2+quadrantMaximum(2,:).^2);
quadrantMaximumNorm = quadrantMaximumNorm';
quadrantMaximumNorm = reshape(quadrantMaximumNorm,16,16);
%imagesc(quadrantMaximumNorm); colormap gray;

end

function meanI = getMeanImage(images)

[h w d] = size(images);
meanI   = sum(images,3)/d;

end